set.seed(123456)
[Needs more text]
In this document we present a work-flow for integration across different omics datasets.
[Note] This is not the final version of the document.
A package with regularized CCA and multiomics DIABLO method is mixOmics. Package igraph is needed for network analysis.
library(mixOmics)
library(igraph)
Package for complex heatmaps.
# BiocManager::install('ComplexHeatmap')
library(ComplexHeatmap)
# https://www.rdocumentation.org/packages/pheatmap/versions/1.0.12/topics/pheatmap
library(pheatmap)
Some original and adapted functions can be found in the file that is silently processed here.
out <- ""
out <- paste(out, knit_child("005-Functions.Rmd", quiet = TRUE))
Usually we use a file management system under the pISA-tree framework. For simplicity, all files ( scripts, data ,… ) are in one directory.
Sample description file (aka phenodata)
pfn <- "../input/phenodata_subset_2023-03-08.txt"
files <- c("../input/data_hormonomics.txt", "../input/data_metabolomics.txt",
"../input/data_qPCR.txt", "../input/data_Phenomics.txt", "../input/data_Proteomics.txt")
For future use and labeling, we need text names of dataset objects.
datanames <- c("hormonomics", "metabolomics", "qPCR", "phenomics", "proteomics")
datanames
## [1] "hormonomics" "metabolomics" "qPCR" "phenomics"
## [5] "proteomics"
Define groups to consider. In a small example, we will pick two groups, based on treatment:
useTreatment <- c("C", "H")
Just to check: Phenodata file name
pfn
## [1] "../input/phenodata_subset_2023-03-08.txt"
Data file names
files
## [1] "../input/data_hormonomics.txt" "../input/data_metabolomics.txt"
## [3] "../input/data_qPCR.txt" "../input/data_Phenomics.txt"
## [5] "../input/data_Proteomics.txt"
It is advisable to first read the phenodata, followed by actual data input. This enables smart selection of samples, based on the sample selection column with the assay name.
pfn
## [1] "../input/phenodata_subset_2023-03-08.txt"
#
phdata <- read.table(pfn, header = TRUE, sep = "\t", stringsAsFactors = FALSE,
row.names = 1)
dim(phdata)
## [1] 56 30
names(phdata)
## [1] "Treatment"
## [2] "Harvest"
## [3] "SamplingDay"
## [4] "DaysOfStressH"
## [5] "DaysOfStressD"
## [6] "DaysOfStressW"
## [7] "DaysRecovery"
## [8] "Replicate"
## [9] "Sample.type"
## [10] "Date"
## [11] "Heat.Drought.Water.Recovery.Days"
## [12] "TreatmentxDatexPlant"
## [13] "TreatmentxSamplingDay"
## [14] "TreatmentxSamplingDayxReplicateNo"
## [15] "Comment"
## [16] "Num_Tubers"
## [17] "Total_Tubers_Weight"
## [18] "FW_SH_Total"
## [19] "DW_SH_Total"
## [20] "Fv.Fm"
## [21] "qL"
## [22] "deltaT"
## [23] "TOP.AREA"
## [24] "COMPACTNESS"
## [25] "WATER.CONSUMPTION"
## [26] "Proteomics"
## [27] "Transcriptomics"
## [28] "Metabolomics"
## [29] "Hormonomics"
## [30] "X_A_300_OverView.R"
pdata <- phdata
Filter out groups
pdata$Treatment
## [1] "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C"
## [14] "C" "C" "C" "D" "D" "D" "D" "D" "D" "D" "D" "H" "H"
## [27] "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H"
## [40] "H" "HD" "HD" "HD" "HD" "HD" "HD" "HD" "HD" "W" "W" "W" "W"
## [53] "W" "W" "W" "W"
pdata <- pdata[pdata$Treatment %in% useTreatment, ]
pdata$Treatment
## [1] "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C"
## [17] "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H" "H"
dim(pdata)
## [1] 32 30
Overview of selected samples:
addmargins(table(pdata$Treatment, pdata$SamplingDay))
##
## 1 7 8 14 Sum
## C 4 4 4 4 16
## H 4 4 4 4 16
## Sum 8 8 8 8 32
.treat <- unique(pdata$Treatment)
.days <- unique(pdata$SamplingDay)
.entry <- 0.5
summary(pdata)
## Treatment Harvest SamplingDay DaysOfStressH
## Length:32 Min. :1.00 Min. : 1.0 Min. : 0.00
## Class :character 1st Qu.:1.75 1st Qu.: 5.5 1st Qu.: 0.00
## Mode :character Median :2.50 Median : 7.5 Median : 0.50
## Mean :2.50 Mean : 7.5 Mean : 3.75
## 3rd Qu.:3.25 3rd Qu.: 9.5 3rd Qu.: 7.25
## Max. :4.00 Max. :14.0 Max. :14.00
## DaysOfStressD DaysOfStressW DaysRecovery Replicate
## Min. :0 Min. :0 Min. :0 Min. :1.00
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:1.75
## Median :0 Median :0 Median :0 Median :2.50
## Mean :0 Mean :0 Mean :0 Mean :2.50
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:3.25
## Max. :0 Max. :0 Max. :0 Max. :4.00
## Sample.type Date
## Length:32 Length:32
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## Heat.Drought.Water.Recovery.Days TreatmentxDatexPlant
## Length:32 Length:32
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## TreatmentxSamplingDay TreatmentxSamplingDayxReplicateNo
## Length:32 Length:32
## Class :character Class :character
## Mode :character Mode :character
##
##
##
## Comment Num_Tubers Total_Tubers_Weight FW_SH_Total
## Length:32 Mode:logical Mode:logical Min. :1
## Class :character NA's:32 NA's:32 1st Qu.:1
## Mode :character Median :1
## Mean :1
## 3rd Qu.:1
## Max. :1
## DW_SH_Total Fv.Fm qL deltaT TOP.AREA
## Min. :1 Min. :1 Min. :1 Min. :1 Min. :1
## 1st Qu.:1 1st Qu.:1 1st Qu.:1 1st Qu.:1 1st Qu.:1
## Median :1 Median :1 Median :1 Median :1 Median :1
## Mean :1 Mean :1 Mean :1 Mean :1 Mean :1
## 3rd Qu.:1 3rd Qu.:1 3rd Qu.:1 3rd Qu.:1 3rd Qu.:1
## Max. :1 Max. :1 Max. :1 Max. :1 Max. :1
## COMPACTNESS WATER.CONSUMPTION Proteomics Transcriptomics
## Min. :1 Min. :1 Min. :1 Min. :1
## 1st Qu.:1 1st Qu.:1 1st Qu.:1 1st Qu.:1
## Median :1 Median :1 Median :1 Median :1
## Mean :1 Mean :1 Mean :1 Mean :1
## 3rd Qu.:1 3rd Qu.:1 3rd Qu.:1 3rd Qu.:1
## Max. :1 Max. :1 Max. :1 Max. :1
## Metabolomics Hormonomics X_A_300_OverView.R
## Min. :1 Min. :1 Min. : 1.0
## 1st Qu.:1 1st Qu.:1 1st Qu.: 5.5
## Median :1 Median :1 Median : 7.5
## Mean :1 Mean :1 Mean : 7.5
## 3rd Qu.:1 3rd Qu.:1 3rd Qu.: 9.5
## Max. :1 Max. :1 Max. :14.0
apply(pdata, 2, function(x) summary(as.factor(x)))
## $Treatment
## C H
## 16 16
##
## $Harvest
## 1 2 3 4
## 8 8 8 8
##
## $SamplingDay
## 1 7 8 14
## 8 8 8 8
##
## $DaysOfStressH
## 0 1 7 8 14
## 16 4 4 4 4
##
## $DaysOfStressD
## 0
## 32
##
## $DaysOfStressW
## 0
## 32
##
## $DaysRecovery
## 0
## 32
##
## $Replicate
## 1 2 3 4
## 8 8 8 8
##
## $Sample.type
## S
## 32
##
## $Date
## 2020-11-04 2020-11-10 2020-11-11 2020-11-17
## 8 8 8 8
##
## $Heat.Drought.Water.Recovery.Days
## 0_0_0_0 1_0_0_0 14_0_0_0 7_0_0_0 8_0_0_0
## 16 4 4 4 4
##
## $TreatmentxDatexPlant
## C_2020-11-04_1 C_2020-11-04_2 C_2020-11-04_3 C_2020-11-04_4
## 1 1 1 1
## C_2020-11-10_1 C_2020-11-10_2 C_2020-11-10_3 C_2020-11-10_4
## 1 1 1 1
## C_2020-11-11_1 C_2020-11-11_2 C_2020-11-11_3 C_2020-11-11_4
## 1 1 1 1
## C_2020-11-17_1 C_2020-11-17_2 C_2020-11-17_3 C_2020-11-17_4
## 1 1 1 1
## H_2020-11-04_1 H_2020-11-04_2 H_2020-11-04_3 H_2020-11-04_4
## 1 1 1 1
## H_2020-11-10_1 H_2020-11-10_2 H_2020-11-10_3 H_2020-11-10_4
## 1 1 1 1
## H_2020-11-11_1 H_2020-11-11_2 H_2020-11-11_3 H_2020-11-11_4
## 1 1 1 1
## H_2020-11-17_1 H_2020-11-17_2 H_2020-11-17_3 H_2020-11-17_4
## 1 1 1 1
##
## $TreatmentxSamplingDay
## C_1 C_14 C_7 C_8 H_1 H_14 H_7 H_8
## 4 4 4 4 4 4 4 4
##
## $TreatmentxSamplingDayxReplicateNo
## C_1_1 C_1_2 C_1_3 C_1_4 C_14_1 C_14_2 C_14_3 C_14_4 C_7_1 C_7_2
## 1 1 1 1 1 1 1 1 1 1
## C_7_3 C_7_4 C_8_1 C_8_2 C_8_3 C_8_4 H_1_1 H_1_2 H_1_3 H_1_4
## 1 1 1 1 1 1 1 1 1 1
## H_14_1 H_14_2 H_14_3 H_14_4 H_7_1 H_7_2 H_7_3 H_7_4 H_8_1 H_8_2
## 1 1 1 1 1 1 1 1 1 1
## H_8_3 H_8_4
## 1 1
##
## $Comment
## only for cor
## 32
##
## $Num_Tubers
## NA's
## 32
##
## $Total_Tubers_Weight
## NA's
## 32
##
## $FW_SH_Total
## 1
## 32
##
## $DW_SH_Total
## 1
## 32
##
## $Fv.Fm
## 1
## 32
##
## $qL
## 1
## 32
##
## $deltaT
## 1
## 32
##
## $TOP.AREA
## 1
## 32
##
## $COMPACTNESS
## 1
## 32
##
## $WATER.CONSUMPTION
## 1
## 32
##
## $Proteomics
## 1
## 32
##
## $Transcriptomics
## 1
## 32
##
## $Metabolomics
## 1
## 32
##
## $Hormonomics
## 1
## 32
##
## $X_A_300_OverView.R
## 1 7 8 14
## 8 8 8 8
For this project we aim to integrate several multi-omics datasets. We have data on hormonomics, metabolomics, and qPCR:
Read all datafiles and assign them to named objects
for (i in 1:length(files)) {
print(datanames[i])
assign(datanames[i], read.table(files[i], header = TRUE, sep = "\t"))
print(head(get(datanames[i])))
}
## [1] "hormonomics"
## SampleID IAA oxIAA IAA.Asp ABA PA DPA
## 1 C_S1_1 37.19565 61.49305 2.224072 36.78497 91.99991 45.63410
## 2 C_S1_2 45.86649 67.84222 1.990994 33.11733 92.61196 62.10651
## 3 C_S1_3 47.60516 52.49759 1.496438 41.66091 93.83119 55.06159
## 4 C_S1_4 37.55562 48.69611 2.240201 39.11995 91.35024 59.82576
## 5 C_S7_1 49.46230 76.79460 1.927162 80.00281 231.68757 178.08837
## 6 C_S7_2 78.91379 93.14059 1.856353 49.76028 166.08573 149.86757
## SA JA JA.Ile X9.10.dhJA X12.OH.JA cisOPDA
## 1 505.2129 2.687303 0.5525212 5.449105 16.13451 353.8608
## 2 519.3793 2.802456 0.5661384 3.700174 23.95458 402.7049
## 3 275.2972 4.761132 0.4273195 2.879761 27.93628 644.9688
## 4 628.0462 4.622662 0.2028412 5.165337 25.72584 750.3532
## 5 1314.9488 6.101443 0.6226329 5.010268 244.12270 1294.9082
## 6 731.2122 3.092070 0.3449961 8.373161 203.74167 873.4217
## [1] "metabolomics"
## SampleID Glukose Fructose Sucrose Starch Asp Glu Asn Ser
## 1 C_S1_1 2.13 2.70 3.449 22.06 1039.5 2514.0 177.8 597.9
## 2 C_S1_2 2.20 2.90 3.382 12.74 844.3 1966.4 167.8 498.1
## 3 C_S1_3 0.82 1.59 2.452 9.98 887.2 2067.8 167.0 441.1
## 4 C_S1_4 2.55 3.01 4.057 15.13 987.6 2347.9 172.3 537.6
## 5 C_S7_1 4.77 3.97 3.533 16.25 793.2 2109.0 276.9 368.5
## 6 C_S7_2 6.30 7.16 6.280 9.57 1391.6 3344.4 498.8 704.2
## Gln Gly His Arg Thr Ala Pro Tyr Val Met Ile Lys
## 1 497.8 137.8 20.7 26.9 225.4 702.5 48.6 25.3 54.1 7.3 48.8 26.0
## 2 408.8 104.6 13.3 29.1 208.5 495.5 57.3 22.1 51.3 6.9 47.0 21.8
## 3 400.1 92.1 17.1 29.5 216.1 514.6 53.5 24.4 52.9 8.6 54.3 26.5
## 4 465.7 117.2 16.8 24.9 252.0 653.1 58.6 22.3 58.4 7.9 52.9 26.7
## 5 265.7 68.7 13.7 38.3 188.6 296.4 85.8 31.2 73.1 4.7 60.4 28.4
## 6 680.9 101.2 19.4 63.9 302.4 664.6 135.7 40.6 96.3 8.6 80.6 33.1
## Leu Phe
## 1 18.8 119.0
## 2 15.7 88.0
## 3 19.5 86.8
## 4 20.1 111.0
## 5 27.4 98.3
## 6 27.1 171.2
## [1] "qPCR"
## SampleID RbohA SnRK2 ACO2 HSP70 PR1b RD29B
## 1 C_S1_1 1.346213 1.497998 0.150379 1.013451 0.324007 0.016995
## 2 C_S1_2 1.496014 1.627412 0.196017 1.067230 0.154344 0.046504
## 3 C_S1_3 1.262812 1.630533 0.440406 0.883357 0.269488 0.016995
## 4 C_S1_4 1.428252 1.530426 0.176572 1.042068 0.255647 0.082387
## 5 C_S7_1 1.404150 1.122600 0.537281 1.028170 0.226714 1.751881
## 6 C_S7_2 0.848707 0.725752 0.150857 0.640827 0.192341 0.113202
## X13.LOX P5CS ERF1 CAT1 CO SWEET SP6A
## 1 0.699873 3.647878 0.560424 0.639920 2.421594 0.816298 0.122374
## 2 0.660413 3.264675 0.637808 0.635107 5.563090 1.873813 0.238811
## 3 0.765764 2.151917 0.585558 0.687279 2.740914 0.930861 0.330151
## 4 0.734167 3.155558 0.663618 0.704880 2.469193 0.933900 0.122374
## 5 1.136201 0.862709 1.634810 0.811006 1.709219 1.494808 3.386457
## 6 0.540666 0.908966 1.239857 0.373886 1.580358 0.853875 0.810631
## M0ZJG3
## 1 1.211115
## 2 1.376397
## 3 1.007479
## 4 0.903495
## 5 2.357865
## 6 1.387190
## [1] "phenomics"
## SampleID FW_SH_Total DW_SH_Total TOP.AREA COMPACTNESS Fv.Fm
## 1 C_S1_1 62.511 2.939 90507.25 0.7167274 0.7065431
## 2 C_S1_2 83.262 4.388 101872.30 0.7121689 0.7133399
## 3 C_S1_3 65.553 3.199 92754.84 0.7108707 0.7137180
## 4 C_S1_4 73.533 3.808 108909.81 0.6882440 0.7064169
## 5 C_S7_1 112.044 7.998 156269.47 0.7066436 0.6953883
## 6 C_S7_2 101.276 6.851 154364.77 0.6910989 0.6940050
## qL deltaT WATER.CONSUMPTION
## 1 1.065883 1.3552322 77
## 2 1.070034 1.2458286 68
## 3 1.089637 0.8838253 78
## 4 1.066492 1.0340176 94
## 5 1.040290 1.2201157 132
## 6 1.031891 1.0551987 133
## [1] "proteomics"
## SampleID VdnDe2_86784 PBdnRY1_7025 VdnPW5_1541 CLCdnPW10_11437
## 1 C_S1_1 0.000413884 0.000469004 0.002308758 NA
## 2 C_S1_2 0.000230027 0.000469191 0.002352451 0.000300619
## 3 C_S1_3 0.000449712 0.000524164 0.001314037 NA
## 4 C_S1_4 0.000255355 0.000520852 0.001661846 0.000222480
## 5 C_S7_1 0.000178758 NA 0.000498582 0.000622978
## 6 C_S7_2 0.000226049 0.000307384 0.000910694 0.000787786
## CLCdnDe2_96512 VdnPW1_81435 CLCdnDe6_4953 CLCdnDe13_5738
## 1 NA 2.82e-05 0.000386853 NA
## 2 NA 2.35e-05 0.000322506 NA
## 3 NA 2.62e-05 0.000520422 NA
## 4 NA NA 0.000397796 NA
## 5 NA NA 0.000222778 NA
## 6 NA NA 0.000563427 NA
## CLCdnPW12_1472 CLCdnPW49_8560 VdnPW4_18801 SdnRY1_10712
## 1 0.000767766 0.000039000 0.000203446 0.000272489
## 2 0.000474884 NA 0.000197874 0.000265026
## 3 0.000530523 0.000072600 0.000189478 0.000338374
## 4 0.000584473 0.000036100 0.000188281 0.000336236
## 5 0.000288815 0.000101007 0.000439347 0.000470757
## 6 0.000432853 NA 0.000407422 0.000347256
## CLCdnDe2_311 CLCdnDe13_63179 CLCdnDe10_24187 TDe_comp110944_c2_seq4
## 1 0.000257351 0.000135326 0.000119484 0.000050600
## 2 0.000178787 0.000037600 0.000033200 NA
## 3 0.000159788 0.000042000 0.000037100 0.000094200
## 4 0.000238167 0.000041700 0.000073700 0.000046800
## 5 0.000222302 0.000116896 NA 0.000131029
## 6 0.000187408 0.000246367 NA 0.000165693
## PBdnRY1_11085 CLCdnPW48_15164 CLCdnPW31_956 VdnDe1_162036
## 1 4.44e-05 0.000263690 0.000334185 0.000146521
## 2 1.85e-05 0.000302265 0.000167159 0.000122150
## 3 4.14e-05 0.000276283 0.000145245 0.000091000
## 4 NA 0.000274538 0.000247419 0.000180799
## 5 NA 0.000427083 0.000288672 0.000253132
## 6 NA 0.000216027 0.000243360 0.000160049
## PBdnRY1_2239 CLCdnPW22_15452 CLCdnPW50_2470 VdnPW3_37011 SdnPW1_46
## 1 0.000165383 0.000161016 NA 0.000123402 1.49e-05
## 2 0.000068900 0.000089500 NA 0.000085700 NA
## 3 0.000057800 0.000100000 NA 0.000057500 NA
## 4 0.000057400 0.000149014 NA 0.000114203 1.38e-05
## 5 0.000107145 NA NA 0.000106596 NA
## 6 0.000045200 0.000058600 NA 0.000112329 3.26e-05
## CLCdnDe7_34340 CLCdnDe11_2286 VdnDe2_53978 CLCdnDe7_6144
## 1 NA NA 0.000323764 0.000176082
## 2 NA NA 0.000337389 0.000366984
## 3 NA NA 0.000452302 0.000245989
## 4 NA NA 0.000524353 0.000244435
## 5 NA 0.000369451 NA 0.000228152
## 6 NA 0.000311459 0.001149384 0.000192340
## CLCdnRY10_6741 PBdnRY1_5389 VdnDe2_29290 VdnDe4_115100
## 1 0.000141861 0.000708110 0.000141561 0.000711383
## 2 0.000059100 0.000641660 0.000094400 0.000773551
## 3 0.000066100 0.000602145 0.000158210 0.000547316
## 4 0.000065600 0.000598341 0.000183412 0.000572483
## 5 NA 0.000877616 0.000073400 0.001041977
## 6 0.000077500 0.000807118 0.000092800 0.001047347
## Sotub11g025210.1.1 VdnDe1_39992 VdnPW4_15664
## 1 0.000228539 0.000389774 NA
## 2 0.000285788 0.000541569 NA
## 3 0.000319272 0.000423515 NA
## 4 0.000211504 0.000360720 NA
## 5 0.000197414 0.000505036 NA
## 6 0.000332853 0.000425761 NA
Check if all 5 objects are created
datanames
## [1] "hormonomics" "metabolomics" "qPCR" "phenomics"
## [5] "proteomics"
datanames %in% ls()
## [1] TRUE TRUE TRUE TRUE TRUE
Datasets for DIABLO need to be collected in a list, with rows corresponding to the same samples. The order of samples from shrinked phenodata will be enforced.
The first component of the list will be a grouping variable, indicating the conditions. We will create reasonable names for groups.
sday <- paste0(0, pdata$SamplingDay)
len <- nchar(sday)
sday <- substr(sday, len - 1, len)
trt <- pdata$Treatment
what <- paste(trt, sday, sep = "")
what
## [1] "C01" "C01" "C01" "C01" "C07" "C07" "C07" "C07" "C08" "C08" "C08"
## [12] "C08" "C14" "C14" "C14" "C14" "H01" "H01" "H01" "H01" "H07" "H07"
## [23] "H07" "H07" "H08" "H08" "H08" "H08" "H14" "H14" "H14" "H14"
X <- list(status = what)
names(X[[1]]) <- rownames(pdata)
X
## $status
## C_S1_1 C_S1_2 C_S1_3 C_S1_4 C_S7_1 C_S7_2 C_S7_3 C_S7_4
## "C01" "C01" "C01" "C01" "C07" "C07" "C07" "C07"
## C_S8_1 C_S8_2 C_S8_3 C_S8_4 C_S14_1 C_S14_2 C_S14_3 C_S14_4
## "C08" "C08" "C08" "C08" "C14" "C14" "C14" "C14"
## H_S1_1 H_S1_2 H_S1_3 H_S1_4 H_S7_1 H_S7_2 H_S7_3 H_S7_4
## "H01" "H01" "H01" "H01" "H07" "H07" "H07" "H07"
## H_S8_1 H_S8_2 H_S8_3 H_S8_4 H_S14_1 H_S14_2 H_S14_3 H_S14_4
## "H08" "H08" "H08" "H08" "H14" "H14" "H14" "H14"
print(addmargins(table(pdata$SamplingDay, what)), zero.print = ".")
## what
## C01 C07 C08 C14 H01 H07 H08 H14 Sum
## 1 4 . . . 4 . . . 8
## 7 . 4 . . . 4 . . 8
## 8 . . 4 . . . 4 . 8
## 14 . . . 4 . . . 4 8
## Sum 4 4 4 4 4 4 4 4 32
print(addmargins(table(pdata$Treatment, what)), zero.print = ".")
## what
## C01 C07 C08 C14 H01 H07 H08 H14 Sum
## C 4 4 4 4 . . . . 16
## H . . . . 4 4 4 4 16
## Sum 4 4 4 4 4 4 4 4 32
Put datasets into the list X and ensure that they all have same order of samples as in phenodata.
datanames
## [1] "hormonomics" "metabolomics" "qPCR" "phenomics"
## [5] "proteomics"
i <- 1
for (i in 1:length(datanames)) {
x <- get(datanames[i])
rownames(x) <- x[, 1]
x <- x[, -1]
X[[i + 1]] <- x[rownames(pdata), ]
names(X)[i + 1] <- datanames[i]
}
str(X)
## List of 6
## $ status : Named chr [1:32] "C01" "C01" "C01" "C01" ...
## ..- attr(*, "names")= chr [1:32] "C_S1_1" "C_S1_2" "C_S1_3" "C_S1_4" ...
## $ hormonomics :'data.frame': 32 obs. of 12 variables:
## ..$ IAA : num [1:32] 37.2 45.9 47.6 37.6 49.5 ...
## ..$ oxIAA : num [1:32] 61.5 67.8 52.5 48.7 76.8 ...
## ..$ IAA.Asp : num [1:32] 2.22 1.99 1.5 2.24 1.93 ...
## ..$ ABA : num [1:32] 36.8 33.1 41.7 39.1 80 ...
## ..$ PA : num [1:32] 92 92.6 93.8 91.4 231.7 ...
## ..$ DPA : num [1:32] 45.6 62.1 55.1 59.8 178.1 ...
## ..$ SA : num [1:32] 505 519 275 628 1315 ...
## ..$ JA : num [1:32] 2.69 2.8 4.76 4.62 6.1 ...
## ..$ JA.Ile : num [1:32] 0.553 0.566 0.427 0.203 0.623 ...
## ..$ X9.10.dhJA: num [1:32] 5.45 3.7 2.88 5.17 5.01 ...
## ..$ X12.OH.JA : num [1:32] 16.1 24 27.9 25.7 244.1 ...
## ..$ cisOPDA : num [1:32] 354 403 645 750 1295 ...
## $ metabolomics:'data.frame': 32 obs. of 22 variables:
## ..$ Glukose : num [1:32] 2.13 2.2 0.82 2.55 4.77 6.3 7.24 3.09 6.34 9.49 ...
## ..$ Fructose: num [1:32] 2.7 2.9 1.59 3.01 3.97 7.16 5.41 4.04 8.52 7.75 ...
## ..$ Sucrose : num [1:32] 3.45 3.38 2.45 4.06 3.53 ...
## ..$ Starch : num [1:32] 22.06 12.74 9.98 15.13 16.25 ...
## ..$ Asp : num [1:32] 1040 844 887 988 793 ...
## ..$ Glu : num [1:32] 2514 1966 2068 2348 2109 ...
## ..$ Asn : num [1:32] 178 168 167 172 277 ...
## ..$ Ser : num [1:32] 598 498 441 538 368 ...
## ..$ Gln : num [1:32] 498 409 400 466 266 ...
## ..$ Gly : num [1:32] 137.8 104.6 92.1 117.2 68.7 ...
## ..$ His : num [1:32] 20.7 13.3 17.1 16.8 13.7 19.4 20.5 17.8 8.8 18.4 ...
## ..$ Arg : num [1:32] 26.9 29.1 29.5 24.9 38.3 63.9 38.6 47.3 54.7 67.2 ...
## ..$ Thr : num [1:32] 225 208 216 252 189 ...
## ..$ Ala : num [1:32] 702 496 515 653 296 ...
## ..$ Pro : num [1:32] 48.6 57.3 53.5 58.6 85.8 ...
## ..$ Tyr : num [1:32] 25.3 22.1 24.4 22.3 31.2 40.6 50.9 37.1 20.2 34.6 ...
## ..$ Val : num [1:32] 54.1 51.3 52.9 58.4 73.1 ...
## ..$ Met : num [1:32] 7.3 6.9 8.6 7.9 4.7 8.6 4.5 6.6 3.8 2.2 ...
## ..$ Ile : num [1:32] 48.8 47 54.3 52.9 60.4 80.6 62.7 63 29.9 48.3 ...
## ..$ Lys : num [1:32] 26 21.8 26.5 26.7 28.4 33.1 29.2 38.4 20.6 41.8 ...
## ..$ Leu : num [1:32] 18.8 15.7 19.5 20.1 27.4 27.1 21.3 26.3 38.5 45.7 ...
## ..$ Phe : num [1:32] 119 88 86.8 111 98.3 ...
## $ qPCR :'data.frame': 32 obs. of 14 variables:
## ..$ RbohA : num [1:32] 1.35 1.5 1.26 1.43 1.4 ...
## ..$ SnRK2 : num [1:32] 1.5 1.63 1.63 1.53 1.12 ...
## ..$ ACO2 : num [1:32] 0.15 0.196 0.44 0.177 0.537 ...
## ..$ HSP70 : num [1:32] 1.013 1.067 0.883 1.042 1.028 ...
## ..$ PR1b : num [1:32] 0.324 0.154 0.269 0.256 0.227 ...
## ..$ RD29B : num [1:32] 0.017 0.0465 0.017 0.0824 1.7519 ...
## ..$ X13.LOX: num [1:32] 0.7 0.66 0.766 0.734 1.136 ...
## ..$ P5CS : num [1:32] 3.648 3.265 2.152 3.156 0.863 ...
## ..$ ERF1 : num [1:32] 0.56 0.638 0.586 0.664 1.635 ...
## ..$ CAT1 : num [1:32] 0.64 0.635 0.687 0.705 0.811 ...
## ..$ CO : num [1:32] 2.42 5.56 2.74 2.47 1.71 ...
## ..$ SWEET : num [1:32] 0.816 1.874 0.931 0.934 1.495 ...
## ..$ SP6A : num [1:32] 0.122 0.239 0.33 0.122 3.386 ...
## ..$ M0ZJG3 : num [1:32] 1.211 1.376 1.007 0.903 2.358 ...
## $ phenomics :'data.frame': 32 obs. of 8 variables:
## ..$ FW_SH_Total : num [1:32] 62.5 83.3 65.6 73.5 112 ...
## ..$ DW_SH_Total : num [1:32] 2.94 4.39 3.2 3.81 8 ...
## ..$ TOP.AREA : num [1:32] 90507 101872 92755 108910 156269 ...
## ..$ COMPACTNESS : num [1:32] 0.717 0.712 0.711 0.688 0.707 ...
## ..$ Fv.Fm : num [1:32] 0.707 0.713 0.714 0.706 0.695 ...
## ..$ qL : num [1:32] 1.07 1.07 1.09 1.07 1.04 ...
## ..$ deltaT : num [1:32] 1.355 1.246 0.884 1.034 1.22 ...
## ..$ WATER.CONSUMPTION: int [1:32] 77 68 78 94 132 133 111 131 144 127 ...
## $ proteomics :'data.frame': 32 obs. of 36 variables:
## ..$ VdnDe2_86784 : num [1:32] 0.000414 0.00023 0.00045 0.000255 0.000179 ...
## ..$ PBdnRY1_7025 : num [1:32] 0.000469 0.000469 0.000524 0.000521 NA ...
## ..$ VdnPW5_1541 : num [1:32] 0.002309 0.002352 0.001314 0.001662 0.000499 ...
## ..$ CLCdnPW10_11437 : num [1:32] NA 0.000301 NA 0.000222 0.000623 ...
## ..$ CLCdnDe2_96512 : num [1:32] NA NA NA NA NA ...
## ..$ VdnPW1_81435 : num [1:32] 2.82e-05 2.35e-05 2.62e-05 NA NA NA NA NA NA NA ...
## ..$ CLCdnDe6_4953 : num [1:32] 0.000387 0.000323 0.00052 0.000398 0.000223 ...
## ..$ CLCdnDe13_5738 : num [1:32] NA NA NA NA NA NA 4.72e-05 NA NA NA ...
## ..$ CLCdnPW12_1472 : num [1:32] 0.000768 0.000475 0.000531 0.000584 0.000289 ...
## ..$ CLCdnPW49_8560 : num [1:32] 3.90e-05 NA 7.26e-05 3.61e-05 1.01e-04 ...
## ..$ VdnPW4_18801 : num [1:32] 0.000203 0.000198 0.000189 0.000188 0.000439 ...
## ..$ SdnRY1_10712 : num [1:32] 0.000272 0.000265 0.000338 0.000336 0.000471 ...
## ..$ CLCdnDe2_311 : num [1:32] 0.000257 0.000179 0.00016 0.000238 0.000222 ...
## ..$ CLCdnDe13_63179 : num [1:32] 1.35e-04 3.76e-05 4.20e-05 4.17e-05 1.17e-04 ...
## ..$ CLCdnDe10_24187 : num [1:32] 1.19e-04 3.32e-05 3.71e-05 7.37e-05 NA ...
## ..$ TDe_comp110944_c2_seq4: num [1:32] 5.06e-05 NA 9.42e-05 4.68e-05 1.31e-04 ...
## ..$ PBdnRY1_11085 : num [1:32] 4.44e-05 1.85e-05 4.14e-05 NA NA NA NA NA NA NA ...
## ..$ CLCdnPW48_15164 : num [1:32] 0.000264 0.000302 0.000276 0.000275 0.000427 ...
## ..$ CLCdnPW31_956 : num [1:32] 0.000334 0.000167 0.000145 0.000247 0.000289 ...
## ..$ VdnDe1_162036 : num [1:32] 0.000147 0.000122 0.000091 0.000181 0.000253 ...
## ..$ PBdnRY1_2239 : num [1:32] 1.65e-04 6.89e-05 5.78e-05 5.74e-05 1.07e-04 ...
## ..$ CLCdnPW22_15452 : num [1:32] 1.61e-04 8.95e-05 1.00e-04 1.49e-04 NA ...
## ..$ CLCdnPW50_2470 : num [1:32] NA NA NA NA NA ...
## ..$ VdnPW3_37011 : num [1:32] 1.23e-04 8.57e-05 5.75e-05 1.14e-04 1.07e-04 ...
## ..$ SdnPW1_46 : num [1:32] 1.49e-05 NA NA 1.38e-05 NA ...
## ..$ CLCdnDe7_34340 : num [1:32] NA NA NA NA NA NA NA NA NA NA ...
## ..$ CLCdnDe11_2286 : num [1:32] NA NA NA NA 0.000369 ...
## ..$ VdnDe2_53978 : num [1:32] 0.000324 0.000337 0.000452 0.000524 NA ...
## ..$ CLCdnDe7_6144 : num [1:32] 0.000176 0.000367 0.000246 0.000244 0.000228 ...
## ..$ CLCdnRY10_6741 : num [1:32] 1.42e-04 5.91e-05 6.61e-05 6.56e-05 NA ...
## ..$ PBdnRY1_5389 : num [1:32] 0.000708 0.000642 0.000602 0.000598 0.000878 ...
## ..$ VdnDe2_29290 : num [1:32] 1.42e-04 9.44e-05 1.58e-04 1.83e-04 7.34e-05 ...
## ..$ VdnDe4_115100 : num [1:32] 0.000711 0.000774 0.000547 0.000572 0.001042 ...
## ..$ Sotub11g025210.1.1 : num [1:32] 0.000229 0.000286 0.000319 0.000212 0.000197 ...
## ..$ VdnDe1_39992 : num [1:32] 0.00039 0.000542 0.000424 0.000361 0.000505 ...
## ..$ VdnPW4_15664 : num [1:32] NA NA NA NA NA NA NA NA 1.66e-05 NA ...
names(X)
## [1] "status" "hormonomics" "metabolomics" "qPCR"
## [5] "phenomics" "proteomics"
Check if sample names in all datasets match.
OK <- TRUE
for (i in 2:length(X)) {
print(ok <- all(names(X[[1]]) == rownames(X[[i]])))
OK <- OK & ok
}
## [1] TRUE
## [1] TRUE
## [1] TRUE
## [1] TRUE
## [1] TRUE
Sample names in datasets match.
Put data into safe object DATA.
DATA <- X
We will also need the names of treatment groups.
groups <- unique(pdata$Treatment)
groups
## [1] "C" "H"
CCDATA <- DATA
names(CCDATA)
## [1] "status" "hormonomics" "metabolomics" "qPCR"
## [5] "phenomics" "proteomics"
out <- ""
out <- paste(out, knit_child("035-DIABLO-2.Rmd", quiet = TRUE))
## Warning in eng_r(options): Failed to tidy R code in chunk 'plotArrow'. Reason:
## Error in base::parse(text = code, keep.source = keep.source) :
## <text>:3:1: unexpected symbol
## 2: title = paste ( groups , collapse = ", " )
## 3: invisible
## ^
cat(out)
Child: 035-DIABLO-2.Rmd ## DIABLO hormonomics, metabolomics, qPCR, phenomics, proteomics ## DIABLO
DIABLO from mixOmics enables
integration of more than two datasets.
Thre datasets are organized as a list of matrices with same samples as rows and variables in columns.
data <- CCDATA[-1]
str(data)
## List of 5
## $ hormonomics :'data.frame': 32 obs. of 12 variables:
## ..$ IAA : num [1:32] 37.2 45.9 47.6 37.6 49.5 ...
## ..$ oxIAA : num [1:32] 61.5 67.8 52.5 48.7 76.8 ...
## ..$ IAA.Asp : num [1:32] 2.22 1.99 1.5 2.24 1.93 ...
## ..$ ABA : num [1:32] 36.8 33.1 41.7 39.1 80 ...
## ..$ PA : num [1:32] 92 92.6 93.8 91.4 231.7 ...
## ..$ DPA : num [1:32] 45.6 62.1 55.1 59.8 178.1 ...
## ..$ SA : num [1:32] 505 519 275 628 1315 ...
## ..$ JA : num [1:32] 2.69 2.8 4.76 4.62 6.1 ...
## ..$ JA.Ile : num [1:32] 0.553 0.566 0.427 0.203 0.623 ...
## ..$ X9.10.dhJA: num [1:32] 5.45 3.7 2.88 5.17 5.01 ...
## ..$ X12.OH.JA : num [1:32] 16.1 24 27.9 25.7 244.1 ...
## ..$ cisOPDA : num [1:32] 354 403 645 750 1295 ...
## $ metabolomics:'data.frame': 32 obs. of 22 variables:
## ..$ Glukose : num [1:32] 2.13 2.2 0.82 2.55 4.77 6.3 7.24 3.09 6.34 9.49 ...
## ..$ Fructose: num [1:32] 2.7 2.9 1.59 3.01 3.97 7.16 5.41 4.04 8.52 7.75 ...
## ..$ Sucrose : num [1:32] 3.45 3.38 2.45 4.06 3.53 ...
## ..$ Starch : num [1:32] 22.06 12.74 9.98 15.13 16.25 ...
## ..$ Asp : num [1:32] 1040 844 887 988 793 ...
## ..$ Glu : num [1:32] 2514 1966 2068 2348 2109 ...
## ..$ Asn : num [1:32] 178 168 167 172 277 ...
## ..$ Ser : num [1:32] 598 498 441 538 368 ...
## ..$ Gln : num [1:32] 498 409 400 466 266 ...
## ..$ Gly : num [1:32] 137.8 104.6 92.1 117.2 68.7 ...
## ..$ His : num [1:32] 20.7 13.3 17.1 16.8 13.7 19.4 20.5 17.8 8.8 18.4 ...
## ..$ Arg : num [1:32] 26.9 29.1 29.5 24.9 38.3 63.9 38.6 47.3 54.7 67.2 ...
## ..$ Thr : num [1:32] 225 208 216 252 189 ...
## ..$ Ala : num [1:32] 702 496 515 653 296 ...
## ..$ Pro : num [1:32] 48.6 57.3 53.5 58.6 85.8 ...
## ..$ Tyr : num [1:32] 25.3 22.1 24.4 22.3 31.2 40.6 50.9 37.1 20.2 34.6 ...
## ..$ Val : num [1:32] 54.1 51.3 52.9 58.4 73.1 ...
## ..$ Met : num [1:32] 7.3 6.9 8.6 7.9 4.7 8.6 4.5 6.6 3.8 2.2 ...
## ..$ Ile : num [1:32] 48.8 47 54.3 52.9 60.4 80.6 62.7 63 29.9 48.3 ...
## ..$ Lys : num [1:32] 26 21.8 26.5 26.7 28.4 33.1 29.2 38.4 20.6 41.8 ...
## ..$ Leu : num [1:32] 18.8 15.7 19.5 20.1 27.4 27.1 21.3 26.3 38.5 45.7 ...
## ..$ Phe : num [1:32] 119 88 86.8 111 98.3 ...
## $ qPCR :'data.frame': 32 obs. of 14 variables:
## ..$ RbohA : num [1:32] 1.35 1.5 1.26 1.43 1.4 ...
## ..$ SnRK2 : num [1:32] 1.5 1.63 1.63 1.53 1.12 ...
## ..$ ACO2 : num [1:32] 0.15 0.196 0.44 0.177 0.537 ...
## ..$ HSP70 : num [1:32] 1.013 1.067 0.883 1.042 1.028 ...
## ..$ PR1b : num [1:32] 0.324 0.154 0.269 0.256 0.227 ...
## ..$ RD29B : num [1:32] 0.017 0.0465 0.017 0.0824 1.7519 ...
## ..$ X13.LOX: num [1:32] 0.7 0.66 0.766 0.734 1.136 ...
## ..$ P5CS : num [1:32] 3.648 3.265 2.152 3.156 0.863 ...
## ..$ ERF1 : num [1:32] 0.56 0.638 0.586 0.664 1.635 ...
## ..$ CAT1 : num [1:32] 0.64 0.635 0.687 0.705 0.811 ...
## ..$ CO : num [1:32] 2.42 5.56 2.74 2.47 1.71 ...
## ..$ SWEET : num [1:32] 0.816 1.874 0.931 0.934 1.495 ...
## ..$ SP6A : num [1:32] 0.122 0.239 0.33 0.122 3.386 ...
## ..$ M0ZJG3 : num [1:32] 1.211 1.376 1.007 0.903 2.358 ...
## $ phenomics :'data.frame': 32 obs. of 8 variables:
## ..$ FW_SH_Total : num [1:32] 62.5 83.3 65.6 73.5 112 ...
## ..$ DW_SH_Total : num [1:32] 2.94 4.39 3.2 3.81 8 ...
## ..$ TOP.AREA : num [1:32] 90507 101872 92755 108910 156269 ...
## ..$ COMPACTNESS : num [1:32] 0.717 0.712 0.711 0.688 0.707 ...
## ..$ Fv.Fm : num [1:32] 0.707 0.713 0.714 0.706 0.695 ...
## ..$ qL : num [1:32] 1.07 1.07 1.09 1.07 1.04 ...
## ..$ deltaT : num [1:32] 1.355 1.246 0.884 1.034 1.22 ...
## ..$ WATER.CONSUMPTION: int [1:32] 77 68 78 94 132 133 111 131 144 127 ...
## $ proteomics :'data.frame': 32 obs. of 36 variables:
## ..$ VdnDe2_86784 : num [1:32] 0.000414 0.00023 0.00045 0.000255 0.000179 ...
## ..$ PBdnRY1_7025 : num [1:32] 0.000469 0.000469 0.000524 0.000521 NA ...
## ..$ VdnPW5_1541 : num [1:32] 0.002309 0.002352 0.001314 0.001662 0.000499 ...
## ..$ CLCdnPW10_11437 : num [1:32] NA 0.000301 NA 0.000222 0.000623 ...
## ..$ CLCdnDe2_96512 : num [1:32] NA NA NA NA NA ...
## ..$ VdnPW1_81435 : num [1:32] 2.82e-05 2.35e-05 2.62e-05 NA NA NA NA NA NA NA ...
## ..$ CLCdnDe6_4953 : num [1:32] 0.000387 0.000323 0.00052 0.000398 0.000223 ...
## ..$ CLCdnDe13_5738 : num [1:32] NA NA NA NA NA NA 4.72e-05 NA NA NA ...
## ..$ CLCdnPW12_1472 : num [1:32] 0.000768 0.000475 0.000531 0.000584 0.000289 ...
## ..$ CLCdnPW49_8560 : num [1:32] 3.90e-05 NA 7.26e-05 3.61e-05 1.01e-04 ...
## ..$ VdnPW4_18801 : num [1:32] 0.000203 0.000198 0.000189 0.000188 0.000439 ...
## ..$ SdnRY1_10712 : num [1:32] 0.000272 0.000265 0.000338 0.000336 0.000471 ...
## ..$ CLCdnDe2_311 : num [1:32] 0.000257 0.000179 0.00016 0.000238 0.000222 ...
## ..$ CLCdnDe13_63179 : num [1:32] 1.35e-04 3.76e-05 4.20e-05 4.17e-05 1.17e-04 ...
## ..$ CLCdnDe10_24187 : num [1:32] 1.19e-04 3.32e-05 3.71e-05 7.37e-05 NA ...
## ..$ TDe_comp110944_c2_seq4: num [1:32] 5.06e-05 NA 9.42e-05 4.68e-05 1.31e-04 ...
## ..$ PBdnRY1_11085 : num [1:32] 4.44e-05 1.85e-05 4.14e-05 NA NA NA NA NA NA NA ...
## ..$ CLCdnPW48_15164 : num [1:32] 0.000264 0.000302 0.000276 0.000275 0.000427 ...
## ..$ CLCdnPW31_956 : num [1:32] 0.000334 0.000167 0.000145 0.000247 0.000289 ...
## ..$ VdnDe1_162036 : num [1:32] 0.000147 0.000122 0.000091 0.000181 0.000253 ...
## ..$ PBdnRY1_2239 : num [1:32] 1.65e-04 6.89e-05 5.78e-05 5.74e-05 1.07e-04 ...
## ..$ CLCdnPW22_15452 : num [1:32] 1.61e-04 8.95e-05 1.00e-04 1.49e-04 NA ...
## ..$ CLCdnPW50_2470 : num [1:32] NA NA NA NA NA ...
## ..$ VdnPW3_37011 : num [1:32] 1.23e-04 8.57e-05 5.75e-05 1.14e-04 1.07e-04 ...
## ..$ SdnPW1_46 : num [1:32] 1.49e-05 NA NA 1.38e-05 NA ...
## ..$ CLCdnDe7_34340 : num [1:32] NA NA NA NA NA NA NA NA NA NA ...
## ..$ CLCdnDe11_2286 : num [1:32] NA NA NA NA 0.000369 ...
## ..$ VdnDe2_53978 : num [1:32] 0.000324 0.000337 0.000452 0.000524 NA ...
## ..$ CLCdnDe7_6144 : num [1:32] 0.000176 0.000367 0.000246 0.000244 0.000228 ...
## ..$ CLCdnRY10_6741 : num [1:32] 1.42e-04 5.91e-05 6.61e-05 6.56e-05 NA ...
## ..$ PBdnRY1_5389 : num [1:32] 0.000708 0.000642 0.000602 0.000598 0.000878 ...
## ..$ VdnDe2_29290 : num [1:32] 1.42e-04 9.44e-05 1.58e-04 1.83e-04 7.34e-05 ...
## ..$ VdnDe4_115100 : num [1:32] 0.000711 0.000774 0.000547 0.000572 0.001042 ...
## ..$ Sotub11g025210.1.1 : num [1:32] 0.000229 0.000286 0.000319 0.000212 0.000197 ...
## ..$ VdnDe1_39992 : num [1:32] 0.00039 0.000542 0.000424 0.000361 0.000505 ...
## ..$ VdnPW4_15664 : num [1:32] NA NA NA NA NA NA NA NA 1.66e-05 NA ...
length(data)
## [1] 5
In addition, outcome, phenotypic state or in our case treatment can also be determined. We have combination of two treatments and four time points.
state <- factor(CCDATA[[1]])
table(state)
## state
## C01 C07 C08 C14 H01 H07 H08 H14
## 4 4 4 4 4 4 4 4
str(state)
## Factor w/ 8 levels "C01","C07","C08",..: 1 1 1 1 2 2 2 2 3 3 ...
## - attr(*, "names")= chr [1:32] "C_S1_1" "C_S1_2" "C_S1_3" "C_S1_4" ...
Note: Additional insights can be found at http://mixomics.org/mixdiablo/diablo-tcga-case-study/.
list.keepX = c(25, 25) # select arbitrary values of features to keep
list.keepY = c(25, 25)
par(mfrow = c(2, 2))
pairs <- combn(1:length(names(data)), 2)
nms <- names(data)
outn <- ""
j <- 1
ncomp <- length(data)
cols <- c("orange1", "brown1", "lightgreen", "lightblue", "pink")[1:length(data)]
pchs <- c(16, 17, 15, 18, 20)[1:length(data)]
j <- 1
cutoff <- 0.5
for (j in 1:ncol(pairs)) {
pair <- pairs[, j]
pair
X <- CCDATA[[pair[1] + 1]]
Y <- CCDATA[[pair[2] + 1]]
list.keepX <- rep(min(ncol(X), 25), ncomp)
list.keepY <- rep(min(ncol(Y), 25), ncomp)
x <- spls(X, Y, ncomp = ncomp, keepX = list.keepX, keepY = list.keepY)
assign(paste0("spls", j), x)
cat("\n", paste(nms[pair]), "\n")
cat("Results in:", paste0("spls", j), "\n")
cat("Correlation between pls variates:\n")
print(round(cor(x$variates$X, x$variates$Y), 5))
#
plotVar(x, cutoff = cutoff, title = paste(nms[pair], collapse = ", "),
legend = c(nms[pair][1], nms[pair][2]), var.names = FALSE, style = "graphics",
pch = pchs[pair], cex = c(2, 2), col = cols[pair])
}
##
## hormonomics metabolomics
## Results in: spls1
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.91126 0.00000 0.00000 0.00000 0.00000
## comp2 0.13029 0.80764 0.00000 0.00000 0.00000
## comp3 0.00165 -0.00098 0.61826 0.00000 0.00000
## comp4 0.01098 0.15104 -0.19116 0.70185 0.00000
## comp5 0.20485 0.21339 0.14759 0.17672 0.74411
##
## hormonomics qPCR
## Results in: spls2
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.69330 0.00000 0.00000 0.00000 0.00000
## comp2 -0.40812 0.63635 0.00000 0.00000 0.00000
## comp3 -0.31053 0.39759 0.54900 0.00000 0.00000
## comp4 -0.06313 0.13686 0.15843 0.56643 0.00000
## comp5 0.01097 -0.10659 -0.09235 -0.44820 0.67564
##
## hormonomics phenomics
## Results in: spls3
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.88537 0.00000 0.00000 0.00000 0.00000
## comp2 -0.09078 0.77298 0.00000 0.00000 0.00000
## comp3 -0.06791 0.15535 0.59832 0.00000 0.00000
## comp4 -0.06005 -0.06360 0.05644 0.73616 0.00000
## comp5 0.17747 -0.35169 -0.36697 0.23815 0.59106
##
## hormonomics proteomics
## Results in: spls4
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.81155 -0.00114 0.00016 0.02427 0.00489
## comp2 0.37042 0.88279 0.00120 -0.01077 0.00641
## comp3 -0.02635 0.12354 0.80974 0.01161 0.00050
## comp4 -0.10965 -0.12524 0.05118 0.70884 -0.00289
## comp5 0.15276 0.19591 -0.12736 -0.30608 0.86428
Circle Correlation Plots for pairwise PLS models on ADAPT data. At most top 25 features for each dimension with correlation above 0.5, are displayed.
##
## metabolomics qPCR
## Results in: spls5
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.74472 0.00000 0.00000 0.00000 0.00000
## comp2 0.38516 0.60146 0.00000 0.00000 0.00000
## comp3 0.27693 0.40129 0.61882 0.00000 0.00000
## comp4 -0.18839 -0.26407 -0.31216 0.42881 0.00000
## comp5 -0.20034 -0.34327 -0.41141 0.50847 0.68746
##
## metabolomics phenomics
## Results in: spls6
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.86102 0.00000 0.00000 0.00000 0.0000
## comp2 0.20665 0.71545 0.00000 0.00000 0.0000
## comp3 0.17444 -0.16699 0.70184 0.00000 0.0000
## comp4 0.22099 0.53084 -0.28532 0.78020 0.0000
## comp5 -0.04168 0.07992 -0.22865 0.02994 0.7715
##
## metabolomics proteomics
## Results in: spls7
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.79399 -0.00814 -0.00062 0.01063 -0.00406
## comp2 0.18909 0.79747 0.00972 0.00378 -0.00054
## comp3 -0.15076 0.06283 0.67921 0.00642 0.00512
## comp4 0.21992 0.30803 0.09922 0.74192 -0.00064
## comp5 -0.05297 0.10507 0.28292 0.05642 0.67331
##
## qPCR phenomics
## Results in: spls8
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.70121 0.00000 0.00000 0.00000 0.00000
## comp2 0.39336 0.70342 0.00000 0.00000 0.00000
## comp3 -0.00676 -0.31594 0.55622 0.00000 0.00000
## comp4 0.31253 0.24391 -0.22860 0.56584 0.00000
## comp5 0.12327 -0.03823 0.17410 0.14778 0.46243
Circle Correlation Plots for pairwise PLS models on ADAPT data. At most top 25 features for each dimension with correlation above 0.5, are displayed.
##
## qPCR proteomics
## Results in: spls9
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.73778 -0.01324 0.00610 0.00301 0.00108
## comp2 0.12129 0.72044 -0.02926 -0.01731 -0.00778
## comp3 0.32872 0.29495 0.65660 -0.00235 -0.00773
## comp4 0.09182 0.22973 0.32202 0.69532 -0.00653
## comp5 -0.12718 -0.03331 -0.01589 0.07374 0.73876
##
## phenomics proteomics
## Results in: spls10
## Correlation between pls variates:
## comp1 comp2 comp3 comp4 comp5
## comp1 0.83424 0.02153 -0.00950 0.00024 -0.01540
## comp2 0.05617 0.80130 -0.01888 0.01789 -0.04537
## comp3 0.09844 0.00513 0.72323 0.00116 0.00177
## comp4 0.25492 -0.14296 0.01300 0.65102 0.00150
## comp5 -0.08345 0.12540 0.09723 -0.24061 0.66352
Circle Correlation Plots for pairwise PLS models on ADAPT data. At most top 25 features for each dimension with correlation above 0.5, are displayed.
Following the suggestion in the source, we will use design matrices with small values. This is supposed to keep low classification error rate.
entry <- .entry
design = matrix(entry, ncol = length(data), nrow = length(data), dimnames = list(names(data),
names(data)))
diag(design) = 0 # set diagonal to 0s
design
## hormonomics metabolomics qPCR phenomics proteomics
## hormonomics 0.0 0.5 0.5 0.5 0.5
## metabolomics 0.5 0.0 0.5 0.5 0.5
## qPCR 0.5 0.5 0.0 0.5 0.5
## phenomics 0.5 0.5 0.5 0.0 0.5
## proteomics 0.5 0.5 0.5 0.5 0.0
With a design in place, the initial DIABLO model can be generated. An arbitrarily high number of components (ncomp = 5) will be used.
# form basic DIABLO model
Y <- state
basic.diablo.model = block.splsda(X = data, Y = Y, ncomp = 5, design = design)
## Design matrix has changed to include Y; each block will be
## linked to Y.
######################################## test eval
knitr::opts_chunk$set(eval = !FALSE)
Details of tuning process can be found in the http://mixomics.org/mixdiablo/diablo-tcga-case-study/.
The process can be computer time consuming and was performed separately.
From previous tuning sessions one can conclude, that the classification rate stays roughly unchanged after two to four components, so we will set the number of components to the number of data sets:
ncomp <- length(data)
ncomp
## [1] 5
We choose the optimal number of variables to select in each data set using the tune.block.splsda() function, for a grid of keepX values for each type of omics. Note that the function has been set to favour a relatively small signature while allowing us to obtain a sufficient number of variables for downstream validation and/or interpretation. See ?tune.block.splsda.
Previous analyses suggest the following list:
$metabolomics
[1] 10 10 5
$hormonomics
[1] 5 5 10
$qPCR
[1] 10 10 5
Number of kept variates is limited with the minimal number of variates in the datasets. We have decided to keep at most 10 variates for each component.
min_variates <- min(sapply(data, ncol), 10)
if (FALSE) keepX <- list(metabolomics = rep(10, ncomp), hormonomics = rep(10,
ncomp), qPCR = rep(10, ncomp))
list.keepX = list()
for (i in 1:length(data)) list.keepX[[i]] <- rep(min_variates, ncomp)
names(list.keepX) <- names(data)
list.keepX
## $hormonomics
## [1] 8 8 8 8 8
##
## $metabolomics
## [1] 8 8 8 8 8
##
## $qPCR
## [1] 8 8 8 8 8
##
## $phenomics
## [1] 8 8 8 8 8
##
## $proteomics
## [1] 8 8 8 8 8
The final DIABLO model is run as:
# set the optimised DIABLO model
final.diablo.model = block.splsda(X = data, Y = as.factor(state), ncomp = ncomp,
keepX = list.keepX, design = design)
## Design matrix has changed to include Y; each block will be
## linked to Y.
The selected variables can be extracted with the function selectVar(), for example in each block, as seen below. Note that the stability of selected variables can be extracted from the output of the perf() function.
# the features selected from components
for (comp in 1:ncomp) {
cat("\nComponent ", comp, ":\n")
for (i in 1:length(data)) {
cat(names(data)[i], "\n")
print(selectVar(final.diablo.model, comp = comp)[[i]]$name)
}
}
##
## Component 1 :
## hormonomics
## [1] "DPA" "SA" "PA" "X12.OH.JA" "X9.10.dhJA"
## [6] "JA.Ile" "ABA" "IAA"
## metabolomics
## [1] "Glukose" "Fructose" "Val" "Tyr" "Ile" "Lys"
## [7] "His" "Gln"
## qPCR
## [1] "X13.LOX" "CAT1" "PR1b" "M0ZJG3" "SP6A" "HSP70"
## [7] "SWEET" "SnRK2"
## phenomics
## [1] "DW_SH_Total" "Fv.Fm" "FW_SH_Total"
## [4] "COMPACTNESS" "deltaT" "WATER.CONSUMPTION"
## [7] "TOP.AREA" "qL"
## proteomics
## [1] "PBdnRY1_2239" "VdnDe1_39992" "CLCdnPW49_8560"
## [4] "PBdnRY1_5389" "CLCdnDe6_4953" "VdnPW3_37011"
## [7] "SdnPW1_46" "PBdnRY1_7025"
##
## Component 2 :
## hormonomics
## [1] "ABA" "cisOPDA" "IAA" "oxIAA" "JA.Ile" "PA"
## [7] "SA" "IAA.Asp"
## metabolomics
## [1] "Starch" "Ser" "Asn" "His" "Arg" "Met"
## [7] "Sucrose" "Gly"
## qPCR
## [1] "SP6A" "RD29B" "SnRK2" "P5CS" "CO" "HSP70" "ACO2" "RbohA"
## phenomics
## [1] "qL" "WATER.CONSUMPTION" "COMPACTNESS"
## [4] "deltaT" "TOP.AREA" "Fv.Fm"
## [7] "FW_SH_Total" "DW_SH_Total"
## proteomics
## [1] "CLCdnPW31_956" "CLCdnPW48_15164"
## [3] "TDe_comp110944_c2_seq4" "VdnPW4_18801"
## [5] "VdnDe2_29290" "VdnDe1_162036"
## [7] "CLCdnRY10_6741" "CLCdnPW22_15452"
##
## Component 3 :
## hormonomics
## [1] "JA.Ile" "oxIAA" "JA" "IAA.Asp" "X9.10.dhJA"
## [6] "DPA" "PA" "SA"
## metabolomics
## [1] "Met" "Gln" "Ala" "Leu" "Glu" "Phe" "Ile" "Arg"
## qPCR
## [1] "ACO2" "PR1b" "RD29B" "CO" "M0ZJG3" "X13.LOX"
## [7] "SP6A" "P5CS"
## phenomics
## [1] "TOP.AREA" "COMPACTNESS" "DW_SH_Total"
## [4] "Fv.Fm" "deltaT" "WATER.CONSUMPTION"
## [7] "FW_SH_Total" "qL"
## proteomics
## [1] "CLCdnPW12_1472" "VdnDe2_53978" "CLCdnDe7_6144"
## [4] "CLCdnDe2_311" "CLCdnPW22_15452" "Sotub11g025210.1.1"
## [7] "VdnDe1_162036" "SdnRY1_10712"
##
## Component 4 :
## hormonomics
## [1] "IAA" "X12.OH.JA" "JA" "cisOPDA" "SA"
## [6] "oxIAA" "ABA" "DPA"
## metabolomics
## [1] "Asp" "Pro" "Sucrose" "Ala" "Glu" "Leu"
## [7] "Phe" "Asn"
## qPCR
## [1] "RbohA" "P5CS" "RD29B" "M0ZJG3" "CO" "SP6A" "HSP70"
## [8] "ERF1"
## phenomics
## [1] "deltaT" "WATER.CONSUMPTION" "FW_SH_Total"
## [4] "COMPACTNESS" "TOP.AREA" "qL"
## [7] "Fv.Fm" "DW_SH_Total"
## proteomics
## [1] "VdnDe2_53978" "VdnDe4_115100"
## [3] "VdnPW5_1541" "VdnPW4_18801"
## [5] "CLCdnPW10_11437" "CLCdnDe13_63179"
## [7] "CLCdnDe10_24187" "TDe_comp110944_c2_seq4"
##
## Component 5 :
## hormonomics
## [1] "cisOPDA" "JA" "oxIAA" "ABA" "X9.10.dhJA"
## [6] "PA" "SA" "X12.OH.JA"
## metabolomics
## [1] "Gly" "Fructose" "Ile" "Glu" "Glukose" "His"
## [7] "Leu" "Tyr"
## qPCR
## [1] "ERF1" "PR1b" "CO" "ACO2" "P5CS" "M0ZJG3"
## [7] "X13.LOX" "RD29B"
## phenomics
## [1] "WATER.CONSUMPTION" "DW_SH_Total" "Fv.Fm"
## [4] "qL" "TOP.AREA" "deltaT"
## [7] "COMPACTNESS" "FW_SH_Total"
## proteomics
## [1] "TDe_comp110944_c2_seq4" "SdnRY1_10712"
## [3] "PBdnRY1_7025" "CLCdnRY10_6741"
## [5] "VdnDe2_86784" "VdnDe1_162036"
## [7] "CLCdnPW10_11437" "VdnPW3_37011"
plotDIABLO() is a diagnostic plot to check whether the correlation between components from each data set has been maximised as specified in the design matrix. We specify which dimension to be assessed with the ncomp argument.
for (comp in 1:ncomp) {
plotDiablo(final.diablo.model, ncomp = comp)
title(paste("Component", comp), adj = 0.1, line = -1, outer = TRUE)
}
The sample plot with the plotIndiv() function projects each sample into the space spanned by the components of each block. Clustering of the samples can be assessed with this plot.
plind <- plotIndiv(final.diablo.model, ind.names = FALSE, legend = TRUE,
title = "DIABLO Sample Plots", guide = "none", ellipse = TRUE)
In the arrow plot below, the start of the arrow indicates the centroid between all data sets for a given sample and the tips of the arrows indicate the location of that sample in each block. Such graphics highlight the agreement between all data sets at the sample level. While somewhat difficult to interpret, even qualitatively, this arrow plot shows proximities of C01 and H01 (both day 1), C07 and C08, and H07 and H08 ( both a day apart). While C samples are in forth quadrant ( D1 < 0, D2 > 0), H samples have ( D1 < 0, D2 < 0) except H14 that is separated on the positive part of Dimension 1.
plotArrow(final.diablo.model, ind.names = FALSE, legend = TRUE,
title = paste(groups,collapse=", ")
)
Several graphical outputs are available to visualise and mine the associations between the selected variables.
The best starting point to evaluate the correlation structure between variables is with the correlation circle plot. A majority of the qPCR variables are positively correlated only with the first component. The hormonomics and metabolomics variables seem to separate along first two dimensions. These first two components correlate highly with the selected variables from the all three dataset. From this, the correlation of each selected feature from all three datasets can be evaluated based on their proximity.
# if(length(data)==3) pick <- 1:3 else pick <- c(4,1:3) cols <-
# c('orange1', 'brown1', 'lightgreen','lightblue')[pick] pchs <-
# c(16, 17, 15, 18)[pick]
cols <- c("orange1", "brown1", "lightgreen", "lightblue", "pink")[1:length(data)]
pchs <- c(16, 17, 15, 18, 20)[1:length(data)]
plotVar(final.diablo.model, var.names = FALSE, style = "graphics", legend = TRUE,
pch = pchs, cex = rep(2, length(data)), col = cols)
The circos plot is exclusive to integrative frameworks and represents the correlations between variables of different types, represented on the side quadrants. It seems that the qPCR variables are almost entirely negatively correlated with the other two dataframes. Just few from metabolomics and hormonomics are positively correlated. Note that these correlations are above a value of 0.7 (cutoff = 0.7). All interpretations made above are only relevant for features with very strong correlations.
Plot variables
# plotVar(res, cutoff=0.5, legend = TRUE, overlap=!FALSE,
# style='graphics') plotVar(res, cutoff=0.5, legend = TRUE,
# overlap=FALSE, style='graphics')
plotVar(final.diablo.model, cutoff = 0.5, legend = FALSE, comp = c(1, 2),
overlap = FALSE, col = cols, cex = rep(4, length(data)))
plotVar(final.diablo.model, cutoff = 0.5, legend = FALSE, comp = c(2, 3),
overlap = FALSE, col = cols, cex = rep(4, length(data)))
circosPlot(final.diablo.model, cutoff = 0.7, line = TRUE, color.blocks = cols,
color.cor = c(3, 2), size.labels = 1, xpd = TRUE)
Another visualisation of the correlations between the different types of variables is the relevance network, which is also built on the similarity matrix (as is the circos plot). Colour represent variable type.
blocks <- combn(length(data), 2)
j <- 1
cutoff <- 0.8
out35a <- ""
for (j in 1:ncol(blocks)) {
out35a <- paste(out35a, knit_child("035a-DIABLO-network.Rmd", quiet = !TRUE))
if (interactive())
readline()
}
cat(out35a)
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-hormonomics-metabolomics-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-hormonomics-metabolomics-8
=======
../output/figs/network-035a-hormonomics-metabolomics-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-hormonomics-qPCR-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-hormonomics-qPCR-8
=======
../output/figs/network-035a-hormonomics-qPCR-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-hormonomics-phenomics-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-hormonomics-phenomics-8
=======
../output/figs/network-035a-hormonomics-phenomics-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-hormonomics-proteomics-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-hormonomics-proteomics-8
=======
../output/figs/network-035a-hormonomics-proteomics-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-metabolomics-qPCR-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-metabolomics-qPCR-8
=======
../output/figs/network-035a-metabolomics-qPCR-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-metabolomics-phenomics-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-metabolomics-phenomics-8
=======
../output/figs/network-035a-metabolomics-phenomics-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-metabolomics-proteomics-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-metabolomics-proteomics-8
=======
../output/figs/network-035a-metabolomics-proteomics-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-qPCR-phenomics-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-qPCR-phenomics-8
=======
../output/figs/network-035a-qPCR-phenomics-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-qPCR-proteomics-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-qPCR-proteomics-8
=======
../output/figs/network-035a-qPCR-proteomics-8
>>>>>>> Stashed changes
nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[,
=======
nfn <- paste0("../output/figs/network-035a-", paste(names(data)[blocks[,
>>>>>>> Stashed changes
j]], collapse = "-"), "-", cutoff * 10)
nfn
## [1] "../reports/figs/network-035a-phenomics-proteomics-8"
png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[,
j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)
Cutoff = 0.8
<<<<<<< Updated upstream
../reports/figs/network-035a-phenomics-proteomics-8
=======
../output/figs/network-035a-phenomics-proteomics-8
>>>>>>> Stashed changes
cutoff <- 0
x <- final.diablo.model
layout.fun <- NULL
label <- paste(.treat, collapse = ", ")
out35b <- ""
out35b <- paste(out35b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
cat(out35b)
ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, cutoff * 10, sep = "-")
set.seed(1234)
<<<<<<< Updated upstream
nw <- network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff,
shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)
## #########
## 035b-multipartite
# nw <- my.network(x , blocks = 1:ndata , color.node = cols , cutoff
# = cutoff , shape.node = 'rectangle' , layout = layout.fun , save =
# 'png' , name.save = nfn )
# Save network layout in ly, used by my.layout function.
if (exists(deparse(substitute(nw)))) ly <- nw$layout else ly <- NULL
cutoff <- 0.8
x <- final.diablo.model
label <- paste(.treat, collapse = ", ")
out35b <- ""
out35b <- paste(out35b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
cat(out35b)
ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, cutoff * 10, sep = "-")
set.seed(1234)
<<<<<<< Updated upstream
nw <- network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff,
shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)
## #########
## 035b-multipartite
# nw <- my.network(x , blocks = 1:ndata , color.node = cols , cutoff
# = cutoff , shape.node = 'rectangle' , layout = layout.fun , save =
# 'png' , name.save = nfn )
The function “plotLoadings” visualises the loading weights of each selected variable on each component and each data set. The colour indicates the class in which the variable has the maximum level of expression “contrib = ‘max’” using the median “method = ‘median’”. Figure below depicts the loading values for each dimension.
for (i in 1:ncomp) plotLoadings(final.diablo.model, comp = i, contrib = "max",
method = "median")
The cimDIABLO() function is a clustered image map specifically implemented to represent the multi-omics molecular signature expression for each sample. From figure below the areas of homogeneous expression levels for a set of samples across a set of features can be determined. For instance, the H14 samples were the only group to show extremely high levels of expression for a specific set of genes and metabolites. This indicates these features are fairly discriminating for this subtype.
cimfn <- "../reports/figs/cim.png"
png(cimfn, res = 600, width = 4000, height = 4000)
cimDiablo(final.diablo.model, size.legend = 0.7)
##
## trimming values to [-3, 3] range for cim visualisation. See 'trim' arg in ?cimDiablo
dev.off()
## pdf
## 2
An AUC plot per block can also be obtained using the function auroc(). The interpretation of this output may not be particularly insightful in relation to the performance evaluation of our methods, but can complement the statistical analysis.
par(mfrow = c(2, 2))
for (i in 1:length(data)) auc.splsda = auroc(final.diablo.model, roc.block = names(data[i]),
roc.comp = 1, print = FALSE)
Save finil DIABLO model for future use in networks.
res <- final.diablo.model
Estimate classification error rate. The error rate should drop by more components used.
# run component number tuning with repeated CV
system.time(perf.diablo <- perf(res, validation = "Mfold", folds = 3, nrepeat = 10))
## user system elapsed
<<<<<<< Updated upstream
## 21.72 0.42 22.58
=======
## 19.69 0.26 19.97
>>>>>>> Stashed changes
plot(perf.diablo) # plot output of tuning
# the features selected to form components
for (comp in 1:ncomp) {
cat("\nComponent ", comp, ":\n")
for (i in 1:length(data)) {
cat(names(data)[i], "\n")
print(selectVar(res, comp = comp)[[i]]$name)
}
}
##
## Component 1 :
## hormonomics
## [1] "DPA" "SA" "PA" "X12.OH.JA" "X9.10.dhJA"
## [6] "JA.Ile" "ABA" "IAA"
## metabolomics
## [1] "Glukose" "Fructose" "Val" "Tyr" "Ile" "Lys"
## [7] "His" "Gln"
## qPCR
## [1] "X13.LOX" "CAT1" "PR1b" "M0ZJG3" "SP6A" "HSP70"
## [7] "SWEET" "SnRK2"
## phenomics
## [1] "DW_SH_Total" "Fv.Fm" "FW_SH_Total"
## [4] "COMPACTNESS" "deltaT" "WATER.CONSUMPTION"
## [7] "TOP.AREA" "qL"
## proteomics
## [1] "PBdnRY1_2239" "VdnDe1_39992" "CLCdnPW49_8560"
## [4] "PBdnRY1_5389" "CLCdnDe6_4953" "VdnPW3_37011"
## [7] "SdnPW1_46" "PBdnRY1_7025"
##
## Component 2 :
## hormonomics
## [1] "ABA" "cisOPDA" "IAA" "oxIAA" "JA.Ile" "PA"
## [7] "SA" "IAA.Asp"
## metabolomics
## [1] "Starch" "Ser" "Asn" "His" "Arg" "Met"
## [7] "Sucrose" "Gly"
## qPCR
## [1] "SP6A" "RD29B" "SnRK2" "P5CS" "CO" "HSP70" "ACO2" "RbohA"
## phenomics
## [1] "qL" "WATER.CONSUMPTION" "COMPACTNESS"
## [4] "deltaT" "TOP.AREA" "Fv.Fm"
## [7] "FW_SH_Total" "DW_SH_Total"
## proteomics
## [1] "CLCdnPW31_956" "CLCdnPW48_15164"
## [3] "TDe_comp110944_c2_seq4" "VdnPW4_18801"
## [5] "VdnDe2_29290" "VdnDe1_162036"
## [7] "CLCdnRY10_6741" "CLCdnPW22_15452"
##
## Component 3 :
## hormonomics
## [1] "JA.Ile" "oxIAA" "JA" "IAA.Asp" "X9.10.dhJA"
## [6] "DPA" "PA" "SA"
## metabolomics
## [1] "Met" "Gln" "Ala" "Leu" "Glu" "Phe" "Ile" "Arg"
## qPCR
## [1] "ACO2" "PR1b" "RD29B" "CO" "M0ZJG3" "X13.LOX"
## [7] "SP6A" "P5CS"
## phenomics
## [1] "TOP.AREA" "COMPACTNESS" "DW_SH_Total"
## [4] "Fv.Fm" "deltaT" "WATER.CONSUMPTION"
## [7] "FW_SH_Total" "qL"
## proteomics
## [1] "CLCdnPW12_1472" "VdnDe2_53978" "CLCdnDe7_6144"
## [4] "CLCdnDe2_311" "CLCdnPW22_15452" "Sotub11g025210.1.1"
## [7] "VdnDe1_162036" "SdnRY1_10712"
##
## Component 4 :
## hormonomics
## [1] "IAA" "X12.OH.JA" "JA" "cisOPDA" "SA"
## [6] "oxIAA" "ABA" "DPA"
## metabolomics
## [1] "Asp" "Pro" "Sucrose" "Ala" "Glu" "Leu"
## [7] "Phe" "Asn"
## qPCR
## [1] "RbohA" "P5CS" "RD29B" "M0ZJG3" "CO" "SP6A" "HSP70"
## [8] "ERF1"
## phenomics
## [1] "deltaT" "WATER.CONSUMPTION" "FW_SH_Total"
## [4] "COMPACTNESS" "TOP.AREA" "qL"
## [7] "Fv.Fm" "DW_SH_Total"
## proteomics
## [1] "VdnDe2_53978" "VdnDe4_115100"
## [3] "VdnPW5_1541" "VdnPW4_18801"
## [5] "CLCdnPW10_11437" "CLCdnDe13_63179"
## [7] "CLCdnDe10_24187" "TDe_comp110944_c2_seq4"
##
## Component 5 :
## hormonomics
## [1] "cisOPDA" "JA" "oxIAA" "ABA" "X9.10.dhJA"
## [6] "PA" "SA" "X12.OH.JA"
## metabolomics
## [1] "Gly" "Fructose" "Ile" "Glu" "Glukose" "His"
## [7] "Leu" "Tyr"
## qPCR
## [1] "ERF1" "PR1b" "CO" "ACO2" "P5CS" "M0ZJG3"
## [7] "X13.LOX" "RD29B"
## phenomics
## [1] "WATER.CONSUMPTION" "DW_SH_Total" "Fv.Fm"
## [4] "qL" "TOP.AREA" "deltaT"
## [7] "COMPACTNESS" "FW_SH_Total"
## proteomics
## [1] "TDe_comp110944_c2_seq4" "SdnRY1_10712"
## [3] "PBdnRY1_7025" "CLCdnRY10_6741"
## [5] "VdnDe2_86784" "VdnDe1_162036"
## [7] "CLCdnPW10_11437" "VdnPW3_37011"
One would like to reduce the number of nodes, especially for proteomics data. One option is to reduce datasets in a way to keep only the variables in the selectVars in original data in . We will keep variables from the first two components.
keptVars <- unique(c(selectVar(res, comp = 1)[[1]]$name, selectVar(res,
comp = 2)[[1]]$name))
which(keptVars %in% selectVar(res, comp = 1)[[1]]$name)
## [1] 1 2 3 4 5 6 7 8
which(keptVars %in% selectVar(res, comp = 2)[[1]]$name)
## [1] 2 3 6 7 8 9 10 11
size.variables <- 1
sim <- circosPlot(final.diablo.model, cutoff = 0.5, line = TRUE, color.blocks = cols,
color.cor = c(3, 2), size.labels = 1, size.variables = size.variables,
xpd = TRUE)
circosPlot(final.diablo.model, cutoff = 0.75, line = TRUE, color.blocks = cols,
color.cor = c(3, 2), size.labels = 1, size.variables = size.variables,
xpd = TRUE)
circosPlot(final.diablo.model, cutoff = 0.9, line = TRUE, color.blocks = cols,
color.cor = c(3, 2), size.labels = 1, size.variables = size.variables,
xpd = TRUE)
We will prepare partial models for each treatment and treatment combination.
netlist <- list()
layout.fun <- NULL
for (i in 1:length(.treat)) {
thisTreat <- .treat[i]
out037 <- ""
out037 <- paste(out037, knit_child("037-network.Rmd", quiet = TRUE))
netlist[[thisTreat]] <- RETURN # complete network
cat(out037)
}
filter <- pdata$Treatment %in% thisTreat
XX1 <- lapply(CCDATA, function(x) if (is.null(dim(x))) x[filter] else x[filter,
])
table(XX1$status)
##
## C01 C07 C08 C14
## 4 4 4 4
res1 <- block.splsda(X = XX1[-1], Y = as.factor(XX1[[1]]), ncomp = ncomp,
keepX = list.keepX, design = design)
## Design matrix has changed to include Y; each block will be
## linked to Y.
cutoff <- 0
x <- res1
layout.fun <- NULL
label <- thisTreat
out23b <- ""
out23b <- paste(out23b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
RETURN <- nw # complete network
cat(out23b)
ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, cutoff * 10, sep = "-")
set.seed(1234)
<<<<<<< Updated upstream
nw <- network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff,
shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)
## #########
## 035b-multipartite
# nw <- my.network(x , blocks = 1:ndata , color.node = cols , cutoff
# = cutoff , shape.node = 'rectangle' , layout = layout.fun , save =
# 'png' , name.save = nfn )
Save network layout for further plots, used by layout function
my.layout.
ly <- nw$layout
cutoff <- 0.7
x <- res1
layout.fun <- my.layout
label <- thisTreat
out23b <- ""
out23b <- paste(out23b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
cat(out23b)
ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, cutoff * 10, sep = "-")
set.seed(1234)
<<<<<<< Updated upstream
nw <- network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff,
shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)
## #########
## 035b-multipartite
# nw <- my.network(x , blocks = 1:ndata , color.node = cols , cutoff
# = cutoff , shape.node = 'rectangle' , layout = layout.fun , save =
# 'png' , name.save = nfn )
filter <- pdata$Treatment %in% thisTreat
XX1 <- lapply(CCDATA, function(x) if (is.null(dim(x))) x[filter] else x[filter,
])
table(XX1$status)
##
## H01 H07 H08 H14
## 4 4 4 4
res1 <- block.splsda(X = XX1[-1], Y = as.factor(XX1[[1]]), ncomp = ncomp,
keepX = list.keepX, design = design)
## Design matrix has changed to include Y; each block will be
## linked to Y.
cutoff <- 0
x <- res1
layout.fun <- NULL
label <- thisTreat
out23b <- ""
out23b <- paste(out23b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
RETURN <- nw # complete network
cat(out23b)
ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, cutoff * 10, sep = "-")
set.seed(1234)
<<<<<<< Updated upstream
nw <- network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff,
shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)
## #########
## 035b-multipartite
# nw <- my.network(x , blocks = 1:ndata , color.node = cols , cutoff
# = cutoff , shape.node = 'rectangle' , layout = layout.fun , save =
# 'png' , name.save = nfn )
Save network layout for further plots, used by layout function
my.layout.
ly <- nw$layout
cutoff <- 0.7
x <- res1
layout.fun <- my.layout
label <- thisTreat
out23b <- ""
out23b <- paste(out23b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
cat(out23b)
ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, cutoff * 10, sep = "-")
set.seed(1234)
<<<<<<< Updated upstream
nw <- network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff,
shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)
## #########
## 035b-multipartite
# nw <- my.network(x , blocks = 1:ndata , color.node = cols , cutoff
# = cutoff , shape.node = 'rectangle' , layout = layout.fun , save =
# 'png' , name.save = nfn )
names(netlist)
## [1] "C" "H"
# Complete network, cutoff = 0, all treatments
datasets <- names(CCDATA[-1])
ndatasets <- length(datasets)
#
nname <- paste("../ouptut/figs/network-", paste(useTreatment, collapse = ""),
sep = "")
N12 <- network(res, cutoff = 0, blocks = 1:ndatasets, shape.node = c("rectangle"),
save = "png", name.save = "network-CH")
#
e <- extractEdges2(N12)
colnames(e)[ncol(e)] <- paste(.treat, collapse = ".")
head(e)
## edge
## ho.IAA_hormonomics_me.Glukose_metabolomics ho.IAA_hormonomics_me.Glukose_metabolomics
## ho.oxIAA_hormonomics_me.Glukose_metabolomics ho.oxIAA_hormonomics_me.Glukose_metabolomics
## ho.IAA.Asp_hormonomics_me.Glukose_metabolomics ho.IAA.Asp_hormonomics_me.Glukose_metabolomics
## ho.ABA_hormonomics_me.Glukose_metabolomics ho.ABA_hormonomics_me.Glukose_metabolomics
## ho.PA_hormonomics_me.Glukose_metabolomics ho.PA_hormonomics_me.Glukose_metabolomics
## ho.DPA_hormonomics_me.Glukose_metabolomics ho.DPA_hormonomics_me.Glukose_metabolomics
## group1
## ho.IAA_hormonomics_me.Glukose_metabolomics hormonomics
## ho.oxIAA_hormonomics_me.Glukose_metabolomics hormonomics
## ho.IAA.Asp_hormonomics_me.Glukose_metabolomics hormonomics
## ho.ABA_hormonomics_me.Glukose_metabolomics hormonomics
## ho.PA_hormonomics_me.Glukose_metabolomics hormonomics
## ho.DPA_hormonomics_me.Glukose_metabolomics hormonomics
## from
## ho.IAA_hormonomics_me.Glukose_metabolomics IAA_hormonomics
## ho.oxIAA_hormonomics_me.Glukose_metabolomics oxIAA_hormonomics
## ho.IAA.Asp_hormonomics_me.Glukose_metabolomics IAA.Asp_hormonomics
## ho.ABA_hormonomics_me.Glukose_metabolomics ABA_hormonomics
## ho.PA_hormonomics_me.Glukose_metabolomics PA_hormonomics
## ho.DPA_hormonomics_me.Glukose_metabolomics DPA_hormonomics
## group2
## ho.IAA_hormonomics_me.Glukose_metabolomics metabolomics
## ho.oxIAA_hormonomics_me.Glukose_metabolomics metabolomics
## ho.IAA.Asp_hormonomics_me.Glukose_metabolomics metabolomics
## ho.ABA_hormonomics_me.Glukose_metabolomics metabolomics
## ho.PA_hormonomics_me.Glukose_metabolomics metabolomics
## ho.DPA_hormonomics_me.Glukose_metabolomics metabolomics
## to
## ho.IAA_hormonomics_me.Glukose_metabolomics Glukose_metabolomics
## ho.oxIAA_hormonomics_me.Glukose_metabolomics Glukose_metabolomics
## ho.IAA.Asp_hormonomics_me.Glukose_metabolomics Glukose_metabolomics
## ho.ABA_hormonomics_me.Glukose_metabolomics Glukose_metabolomics
## ho.PA_hormonomics_me.Glukose_metabolomics Glukose_metabolomics
## ho.DPA_hormonomics_me.Glukose_metabolomics Glukose_metabolomics
## C.H
## ho.IAA_hormonomics_me.Glukose_metabolomics -0.0448342
## ho.oxIAA_hormonomics_me.Glukose_metabolomics 0.4602270
## ho.IAA.Asp_hormonomics_me.Glukose_metabolomics -0.1825575
## ho.ABA_hormonomics_me.Glukose_metabolomics 0.7655279
## ho.PA_hormonomics_me.Glukose_metabolomics 0.8893326
## ho.DPA_hormonomics_me.Glukose_metabolomics 0.7218137
tail(e)
## edge
## ph.TOP.AREA_phenomics_pr.VdnDe1_39992_proteomics ph.TOP.AREA_phenomics_pr.VdnDe1_39992_proteomics
## ph.COMPACTNESS_phenomics_pr.VdnDe1_39992_proteomics ph.COMPACTNESS_phenomics_pr.VdnDe1_39992_proteomics
## ph.Fv.Fm_phenomics_pr.VdnDe1_39992_proteomics ph.Fv.Fm_phenomics_pr.VdnDe1_39992_proteomics
## ph.qL_phenomics_pr.VdnDe1_39992_proteomics ph.qL_phenomics_pr.VdnDe1_39992_proteomics
## ph.deltaT_phenomics_pr.VdnDe1_39992_proteomics ph.deltaT_phenomics_pr.VdnDe1_39992_proteomics
## ph.WATER.CONSUMPTION_phenomics_pr.VdnDe1_39992_proteomics ph.WATER.CONSUMPTION_phenomics_pr.VdnDe1_39992_proteomics
## group1
## ph.TOP.AREA_phenomics_pr.VdnDe1_39992_proteomics phenomics
## ph.COMPACTNESS_phenomics_pr.VdnDe1_39992_proteomics phenomics
## ph.Fv.Fm_phenomics_pr.VdnDe1_39992_proteomics phenomics
## ph.qL_phenomics_pr.VdnDe1_39992_proteomics phenomics
## ph.deltaT_phenomics_pr.VdnDe1_39992_proteomics phenomics
## ph.WATER.CONSUMPTION_phenomics_pr.VdnDe1_39992_proteomics phenomics
## from
## ph.TOP.AREA_phenomics_pr.VdnDe1_39992_proteomics TOP.AREA_phenomics
## ph.COMPACTNESS_phenomics_pr.VdnDe1_39992_proteomics COMPACTNESS_phenomics
## ph.Fv.Fm_phenomics_pr.VdnDe1_39992_proteomics Fv.Fm_phenomics
## ph.qL_phenomics_pr.VdnDe1_39992_proteomics qL_phenomics
## ph.deltaT_phenomics_pr.VdnDe1_39992_proteomics deltaT_phenomics
## ph.WATER.CONSUMPTION_phenomics_pr.VdnDe1_39992_proteomics WATER.CONSUMPTION_phenomics
## group2
## ph.TOP.AREA_phenomics_pr.VdnDe1_39992_proteomics proteomics
## ph.COMPACTNESS_phenomics_pr.VdnDe1_39992_proteomics proteomics
## ph.Fv.Fm_phenomics_pr.VdnDe1_39992_proteomics proteomics
## ph.qL_phenomics_pr.VdnDe1_39992_proteomics proteomics
## ph.deltaT_phenomics_pr.VdnDe1_39992_proteomics proteomics
## ph.WATER.CONSUMPTION_phenomics_pr.VdnDe1_39992_proteomics proteomics
## to
## ph.TOP.AREA_phenomics_pr.VdnDe1_39992_proteomics VdnDe1_39992_proteomics
## ph.COMPACTNESS_phenomics_pr.VdnDe1_39992_proteomics VdnDe1_39992_proteomics
## ph.Fv.Fm_phenomics_pr.VdnDe1_39992_proteomics VdnDe1_39992_proteomics
## ph.qL_phenomics_pr.VdnDe1_39992_proteomics VdnDe1_39992_proteomics
## ph.deltaT_phenomics_pr.VdnDe1_39992_proteomics VdnDe1_39992_proteomics
## ph.WATER.CONSUMPTION_phenomics_pr.VdnDe1_39992_proteomics VdnDe1_39992_proteomics
## C.H
## ph.TOP.AREA_phenomics_pr.VdnDe1_39992_proteomics -0.40647170
## ph.COMPACTNESS_phenomics_pr.VdnDe1_39992_proteomics -0.42836434
## ph.Fv.Fm_phenomics_pr.VdnDe1_39992_proteomics 0.68250477
## ph.qL_phenomics_pr.VdnDe1_39992_proteomics 0.03964368
## ph.deltaT_phenomics_pr.VdnDe1_39992_proteomics -0.48231592
## ph.WATER.CONSUMPTION_phenomics_pr.VdnDe1_39992_proteomics -0.26855605
dim(e)
## [1] 2630 6
Save networks file for combined and single treatments, stored in
netlist. Combine edges from all networks as columns in a
data.frame.
for (i in 1:length(netlist)) {
e1 <- extractEdges2(netlist[[i]])
colnames(e1)[ncol(e1)] <- names(netlist)[i]
e <- merge(e, e1, sort = FALSE, all = TRUE)
}
str(e)
## 'data.frame': 2630 obs. of 8 variables:
## $ edge : chr "ph.FW_SH_Total_phenomics_pr.VdnDe2_86784_proteomics" "ph.DW_SH_Total_phenomics_pr.VdnDe2_86784_proteomics" "ph.TOP.AREA_phenomics_pr.VdnDe2_86784_proteomics" "ph.COMPACTNESS_phenomics_pr.VdnDe2_86784_proteomics" ...
## $ group1: chr "phenomics" "phenomics" "phenomics" "phenomics" ...
## $ from : chr "FW_SH_Total_phenomics" "DW_SH_Total_phenomics" "TOP.AREA_phenomics" "COMPACTNESS_phenomics" ...
## $ group2: chr "proteomics" "proteomics" "proteomics" "proteomics" ...
## $ to : chr "VdnDe2_86784_proteomics" "VdnDe2_86784_proteomics" "VdnDe2_86784_proteomics" "VdnDe2_86784_proteomics" ...
## $ C.H : num -0.0999 -0.2028 0.1263 -0.4047 0.2552 ...
## $ C : num 0 0 0 0 0 0 0 0 0 0 ...
## $ H : num 0 0 0 0 0 0 0 0 0 0 ...
Compose file name and necessary information for network export file
netfn <- paste0("../output/network-", paste(.treat, collapse = "_"), "-",
paste(datasets, collapse = "_"), ".txt")
label0 <- paste(paste(.treat, collapse = ", "), "|", paste(datasets, collapse = ", "),
"; cutoff =", 0)
title <- label0
sets <- 1:length(DATA)
suffix <- paste0(substr(names(DATA), 1, 2)[sets[-1]], collapse = "-")
# e <- data.frame(x=1:10,y=1:10) my.write.table(e,
# file='network.txt',meta=FALSE)
write.table(e, file = netfn, na = "0", sep = "\t", col.names = NA, row.names = TRUE)
Network edges are exported into file
netfn
## [1] "../output/network-C_H-hormonomics_metabolomics_qPCR_phenomics_proteomics.txt"
Table with edges for networks based on combined treatments (C, H) and single treatments (C) and (H) is exported as a text file. This table can be used for inspection and filtering out edges based on selected cutoff. Missing edges are labeled as weight 0. This enables numeric filtration in other visualization or analysis files (e.g. Excel).
Windows 10 x64 (build 19045)
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64 Running under: Windows 10 x64 (build
19045)
Matrix products: default
locale: [1] LC_COLLATE=English_United Kingdom.utf8 [2]
LC_CTYPE=English_United Kingdom.utf8
[3] LC_MONETARY=English_United Kingdom.utf8 [4] LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.utf8
time zone: Europe/Ljubljana tzcode source: internal
attached base packages: [1] grid stats graphics grDevices utils datasets [7] methods base
other attached packages: [1] pheatmap_1.0.12 ComplexHeatmap_2.22.0
igraph_2.1.1
[4] mixOmics_6.30.0 ggplot2_3.5.1 lattice_0.22-6
[7] MASS_7.3-61 knitr_1.48
loaded via a namespace (and not attached): [1] shape_1.4.6.1
circlize_0.4.16 gtable_0.3.6
[4] ellipse_0.5.0 rjson_0.2.23 xfun_0.49
[7] bslib_0.8.0 GlobalOptions_0.1.2 ggrepel_0.9.6
[10] vctrs_0.6.5 tools_4.4.1 generics_0.1.3
[13] stats4_4.4.1 parallel_4.4.1 tibble_3.2.1
[16] fansi_1.0.6 highr_0.11 cluster_2.1.6
[19] rARPACK_0.11-0 pkgconfig_2.0.3 Matrix_1.7-1
[22] RColorBrewer_1.1-3 S4Vectors_0.44.0 lifecycle_1.0.4
[25] farver_2.1.2 compiler_4.4.1 stringr_1.5.1
[28] munsell_0.5.1 codetools_0.2-20 clue_0.3-65
[31] htmltools_0.5.8.1 sass_0.4.9 yaml_2.3.10
[34] pillar_1.9.0 crayon_1.5.3 jquerylib_0.1.4
[37] tidyr_1.3.1 BiocParallel_1.40.0 cachem_1.1.0
[40] iterators_1.0.14 foreach_1.5.2 RSpectra_0.16-2
[43] tidyselect_1.2.1 digest_0.6.37 stringi_1.8.4
[46] dplyr_1.1.4 reshape2_1.4.4 purrr_1.0.2
[49] labeling_0.4.3 fastmap_1.2.0 colorspace_2.1-1
[52] cli_3.6.3 magrittr_2.0.3 utf8_1.2.4
[55] corpcor_1.6.10 withr_3.0.2 scales_1.3.0
[58] rmarkdown_2.29 matrixStats_1.4.1 gridExtra_2.3
[61] png_0.1-8 GetoptLong_1.0.5 evaluate_1.0.1
[64] IRanges_2.40.0 doParallel_1.0.17 rlang_1.1.4
[67] Rcpp_1.0.13-1 glue_1.8.0 formatR_1.14
[70] BiocGenerics_0.52.0 rstudioapi_0.17.1 jsonlite_1.8.9
[73] R6_2.5.1 plyr_1.8.9